library(ggplot2)
library(tidyverse)
library(Seurat)
library(miscTools)
library(viridis)
library(venn)

Read in mutation table

head(mutation_table)

head(mutation_table)
length(intersect(rownames(mutation_table), rownames(tiss_subset_tumor@meta.data)))
[1] 3620
#drop X.1 and X.2 to clean mutation table
mutation_table$X.1 <- NULL
mutation_table$X.2 <- NULL

Normalize the number of mutations by the raw # of transcripts found

# correct colnames
colnames(mutation_table) <- gsub(pattern = "\\.", replacement = "-", x = colnames(mutation_table))
# Subset table to tumor cells alone and 
cells.to.use.tumor <- tiss_subset_tumor@meta.data$cell_id
mutation_table <- mutation_table[cells.to.use.tumor, ]
# save the names of the mutated genes
genes.to.use.mutations <- colnames(mutation_table)
length(genes.to.use.mutations)
[1] 14397
# set raw counts equal to the genes in the mutation table and subset to only tumor cells
raw_counts <- as.data.frame(as.matrix(tiss_subset_tumor@raw.data))[genes.to.use.mutations, cells.to.use.tumor]
# transpose the raw counts so that colnames are row ids and rownames are cell ids (same as the mutation table)
raw_counts1 <- as.data.frame(t(raw_counts))
dim(raw_counts1)
[1]  3620 14397
# normalize the mutation table by the raw counts
normalized_mutation_table <- mutation_table/raw_counts1
# Change NANs, NAs or Infs to 0
is.nan.data.frame <- function(x)
do.call(cbind, lapply(x, is.nan))
is.infinite.data.frame <- function(x)
do.call(cbind, lapply(x, is.infinite))
normalized_mutation_table[is.nan(normalized_mutation_table)] <- 0
normalized_mutation_table[is.na(normalized_mutation_table)] <- 0
normalized_mutation_table[is.infinite(normalized_mutation_table)] <- 0

Total Mutations Analysis

# Set Table to Report Only the Total Number of Mutations
total_mutations_table <- as.data.frame(rowSums(normalized_mutation_table))
head(total_mutations_table)
total_mutations_table$cell_name <- rownames(total_mutations_table)
colnames(total_mutations_table) <- c('total_mutations_score','cell_id')
head(total_mutations_table)

Add Metadata to normalized mutation table

Total Mutations detected by Group

pr_mutation_table <- filter(total_mutations_table_wmetadata, analysis == 'grouped_pr')
pd_mutation_table <- filter(total_mutations_table_wmetadata, analysis == 'grouped_pd')
naive_mutation_table <- filter(total_mutations_table_wmetadata, analysis == 'naive')
# Driver Gene
total_mutations_table_driver <- filter(total_mutations_table_wmetadata, driver_gene == 'EGFR' | driver_gene == 'ALK')

By Group

pdf(file = paste(dir, "plot_out/NI08/NI08_total_muts_by_group.pdf", sep = ""))
ggplot(total_mutations_table_wmetadata, aes(x = analysis, y = total_mutations_score, fill = analysis)) + geom_boxplot(outlier.shape = NA) + 
    geom_jitter(alpha = 1/8, aes(colour = analysis)) + guides(colour = FALSE, fill = FALSE) + xlab("Group") + ylab("Mutation Score") + 
    ggtitle("Distribution of Mutation Scores per Cell") + 
    geom_signif(comparisons = list(c("grouped_pd", "grouped_pr")), map_signif_level=TRUE, y_position = 45, annotations = '< 2e-16') + 
    geom_signif(comparisons = list(c("grouped_pd", "naive")), map_signif_level=TRUE, y_position = 50, annotations = '< 2e-16') + 
    geom_signif(comparisons = list(c("grouped_pr", "naive")), map_signif_level=TRUE, y_position = 40, annotations = '0.023') +
    scale_x_discrete(labels = c("naive" = "TN", "grouped_pr" = "PER", "grouped_pd" = "PD"), limits=c("naive","grouped_pr","grouped_pd"))
dev.off()
null device 
          1 
pdf(file = paste(dir, "plot_out/NI08/NI08_total_muts_by_group_ridgeplot.pdf", sep = ""))
ggplot(total_mutations_table_wmetadata, aes(x = total_mutations_score, y = analysis, fill = analysis)) + geom_density_ridges(alpha=.7) + 
    ylab("Group") + xlab("Mutation Score") + scale_y_discrete(labels = c("PD", "PER", "TN")) + 
    theme(legend.position = "none") + theme(axis.title.y = element_blank())
dev.off()
null device 
          1 
pdf(file = paste(dir, "plot_out/NI08/NI08_total_muts_by_group_and_patient.pdf", sep = ""))
ggplot(total_mutations_table_wmetadata, aes(x = analysis, y = total_mutations_score, color = patient_id)) +
    geom_point(position = position_jitterdodge(dodge.width=0.9)) +
    geom_boxplot(outlier.colour = NA, position = position_dodge(width=0.9)) +
    xlab("Group") +
    ylab("Mutation Score") +
    ggtitle("Distribution of Mutation Scores per Cell") +
    guides(colour=guide_legend(title="Response Group")) +
    guides(fill = FALSE, color = FALSE) +
    scale_x_discrete(labels = c("naive" = "TN", "grouped_pr" = "PER", "grouped_pd" = "PD"), limits=c("naive","grouped_pr","grouped_pd"))
dev.off()
null device 
          1 
pdf(file = paste(dir, "plot_out/NI08/NI08_total_muts_by_primvsmet.pdf", sep = ""))
ggplot(data=subset(total_mutations_table_wmetadata, !is.na(primary_or_metastaic)), aes(x = primary_or_metastaic , y = total_mutations_score, color = analysis)) +
      geom_point(position = position_jitterdodge(dodge.width=0.9)) +
      geom_boxplot(outlier.colour = NA, position = position_dodge(width=0.9)) +
      xlab("Group") +
      ylab("Mutation Score") +
      ggtitle("Distribution of Mutation Scores per Cell") +
      guides(colour=guide_legend(title="Response Group")) +
      guides(fill = FALSE)
dev.off()
null device 
          1 
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoU2V1cmF0KQpsaWJyYXJ5KG1pc2NUb29scykKbGlicmFyeSh2aXJpZGlzKQpsaWJyYXJ5KHZlbm4pCmBgYAoKUmVhZCBpbiBtdXRhdGlvbiB0YWJsZQpgYGB7cn0KIyBybShsaXN0PWxzKCkpCmRpciA8LSAiL215Vm9sdW1lL3NjZWxsX2x1bmdfYWRlbm9jYXJjaW5vbWEvIgpsb2FkKGZpbGU9cGFzdGUoZGlyLCJEYXRhX2lucHV0L05JMDRfdHVtb3Jfc2V1cmF0X29iamVjdC5SRGF0YSIsIHNlcD0iIikpCm11dGF0aW9uX3RhYmxlIDwtIHJlYWQuY3N2KGZpbGUgPSBwYXN0ZShkaXIsICJEYXRhX2lucHV0L2dlbmVDZWxsTXV0YXRpb25Db3VudHNfMi4yNy4xOS5jc3YiLCBzZXAgPSAiIiksIGhlYWRlciA9IFQsIHJvdy5uYW1lcyA9IDEpCmhlYWQobXV0YXRpb25fdGFibGUpCmxlbmd0aChpbnRlcnNlY3Qocm93bmFtZXMobXV0YXRpb25fdGFibGUpLCByb3duYW1lcyh0aXNzX3N1YnNldF90dW1vckBtZXRhLmRhdGEpKSkKYGBgCgpgYGB7cn0KI2Ryb3AgWC4xIGFuZCBYLjIgdG8gY2xlYW4gbXV0YXRpb24gdGFibGUKbXV0YXRpb25fdGFibGUkWC4xIDwtIE5VTEwKbXV0YXRpb25fdGFibGUkWC4yIDwtIE5VTEwKYGBgCgpOb3JtYWxpemUgdGhlIG51bWJlciBvZiBtdXRhdGlvbnMgYnkgdGhlIHJhdyAjIG9mIHRyYW5zY3JpcHRzIGZvdW5kCmBgYHtyfQojIGNvcnJlY3QgY29sbmFtZXMKY29sbmFtZXMobXV0YXRpb25fdGFibGUpIDwtIGdzdWIocGF0dGVybiA9ICJcXC4iLCByZXBsYWNlbWVudCA9ICItIiwgeCA9IGNvbG5hbWVzKG11dGF0aW9uX3RhYmxlKSkKIyBTdWJzZXQgdGFibGUgdG8gdHVtb3IgY2VsbHMgYWxvbmUgYW5kIApjZWxscy50by51c2UudHVtb3IgPC0gdGlzc19zdWJzZXRfdHVtb3JAbWV0YS5kYXRhJGNlbGxfaWQKbXV0YXRpb25fdGFibGUgPC0gbXV0YXRpb25fdGFibGVbY2VsbHMudG8udXNlLnR1bW9yLCBdCiMgc2F2ZSB0aGUgbmFtZXMgb2YgdGhlIG11dGF0ZWQgZ2VuZXMKZ2VuZXMudG8udXNlLm11dGF0aW9ucyA8LSBjb2xuYW1lcyhtdXRhdGlvbl90YWJsZSkKbGVuZ3RoKGdlbmVzLnRvLnVzZS5tdXRhdGlvbnMpCiMgc2V0IHJhdyBjb3VudHMgZXF1YWwgdG8gdGhlIGdlbmVzIGluIHRoZSBtdXRhdGlvbiB0YWJsZSBhbmQgc3Vic2V0IHRvIG9ubHkgdHVtb3IgY2VsbHMKcmF3X2NvdW50cyA8LSBhcy5kYXRhLmZyYW1lKGFzLm1hdHJpeCh0aXNzX3N1YnNldF90dW1vckByYXcuZGF0YSkpW2dlbmVzLnRvLnVzZS5tdXRhdGlvbnMsIGNlbGxzLnRvLnVzZS50dW1vcl0KIyB0cmFuc3Bvc2UgdGhlIHJhdyBjb3VudHMgc28gdGhhdCBjb2xuYW1lcyBhcmUgcm93IGlkcyBhbmQgcm93bmFtZXMgYXJlIGNlbGwgaWRzIChzYW1lIGFzIHRoZSBtdXRhdGlvbiB0YWJsZSkKcmF3X2NvdW50czEgPC0gYXMuZGF0YS5mcmFtZSh0KHJhd19jb3VudHMpKQpkaW0ocmF3X2NvdW50czEpCiMgbm9ybWFsaXplIHRoZSBtdXRhdGlvbiB0YWJsZSBieSB0aGUgcmF3IGNvdW50cwpub3JtYWxpemVkX211dGF0aW9uX3RhYmxlIDwtIG11dGF0aW9uX3RhYmxlL3Jhd19jb3VudHMxCiMgQ2hhbmdlIE5BTnMsIE5BcyBvciBJbmZzIHRvIDAKaXMubmFuLmRhdGEuZnJhbWUgPC0gZnVuY3Rpb24oeCkKZG8uY2FsbChjYmluZCwgbGFwcGx5KHgsIGlzLm5hbikpCmlzLmluZmluaXRlLmRhdGEuZnJhbWUgPC0gZnVuY3Rpb24oeCkKZG8uY2FsbChjYmluZCwgbGFwcGx5KHgsIGlzLmluZmluaXRlKSkKbm9ybWFsaXplZF9tdXRhdGlvbl90YWJsZVtpcy5uYW4obm9ybWFsaXplZF9tdXRhdGlvbl90YWJsZSldIDwtIDAKbm9ybWFsaXplZF9tdXRhdGlvbl90YWJsZVtpcy5uYShub3JtYWxpemVkX211dGF0aW9uX3RhYmxlKV0gPC0gMApub3JtYWxpemVkX211dGF0aW9uX3RhYmxlW2lzLmluZmluaXRlKG5vcm1hbGl6ZWRfbXV0YXRpb25fdGFibGUpXSA8LSAwCmBgYAoKVG90YWwgTXV0YXRpb25zIEFuYWx5c2lzCmBgYHtyfQojIFNldCBUYWJsZSB0byBSZXBvcnQgT25seSB0aGUgVG90YWwgTnVtYmVyIG9mIE11dGF0aW9ucwp0b3RhbF9tdXRhdGlvbnNfdGFibGUgPC0gYXMuZGF0YS5mcmFtZShyb3dTdW1zKG5vcm1hbGl6ZWRfbXV0YXRpb25fdGFibGUpKQpoZWFkKHRvdGFsX211dGF0aW9uc190YWJsZSkKdG90YWxfbXV0YXRpb25zX3RhYmxlJGNlbGxfbmFtZSA8LSByb3duYW1lcyh0b3RhbF9tdXRhdGlvbnNfdGFibGUpCmNvbG5hbWVzKHRvdGFsX211dGF0aW9uc190YWJsZSkgPC0gYygndG90YWxfbXV0YXRpb25zX3Njb3JlJywnY2VsbF9pZCcpCmhlYWQodG90YWxfbXV0YXRpb25zX3RhYmxlKQpgYGAKCkFkZCBNZXRhZGF0YSB0byBub3JtYWxpemVkIG11dGF0aW9uIHRhYmxlCmBgYHtyfQp0b3RhbF9tdXRhdGlvbnNfdGFibGUkY2VsbF9pZCA8LSByb3duYW1lcyh0b3RhbF9tdXRhdGlvbnNfdGFibGUpCnRvdGFsX211dGF0aW9uc190YWJsZV93bWV0YWRhdGEgPC0gbGVmdF9qb2luKHRvdGFsX211dGF0aW9uc190YWJsZSwgdGlzc19zdWJzZXRfdHVtb3JAbWV0YS5kYXRhLCBieSA9ICdjZWxsX2lkJykKcm93bmFtZXModG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSkgPC0gdG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSRjZWxsX2lkCmBgYAoKVG90YWwgTXV0YXRpb25zIGRldGVjdGVkIGJ5IEdyb3VwCmBgYHtyfQpwcl9tdXRhdGlvbl90YWJsZSA8LSBmaWx0ZXIodG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSwgYW5hbHlzaXMgPT0gJ2dyb3VwZWRfcHInKQpwZF9tdXRhdGlvbl90YWJsZSA8LSBmaWx0ZXIodG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSwgYW5hbHlzaXMgPT0gJ2dyb3VwZWRfcGQnKQpuYWl2ZV9tdXRhdGlvbl90YWJsZSA8LSBmaWx0ZXIodG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSwgYW5hbHlzaXMgPT0gJ25haXZlJykKCiMgRHJpdmVyIEdlbmUKdG90YWxfbXV0YXRpb25zX3RhYmxlX2RyaXZlciA8LSBmaWx0ZXIodG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSwgZHJpdmVyX2dlbmUgPT0gJ0VHRlInIHwgZHJpdmVyX2dlbmUgPT0gJ0FMSycpCmBgYAoKQnkgR3JvdXAKYGBge3J9CnRvdGFsX210X3Njb3JlIDwtIHBhaXJ3aXNlLndpbGNveC50ZXN0KHggPSB0b3RhbF9tdXRhdGlvbnNfdGFibGVfd21ldGFkYXRhJHRvdGFsX211dGF0aW9uc19zY29yZSwgZyA9IHRvdGFsX211dGF0aW9uc190YWJsZV93bWV0YWRhdGEkYW5hbHlzaXMpCgpwZGYoZmlsZSA9IHBhc3RlKGRpciwgInBsb3Rfb3V0L05JMDgvTkkwOF90b3RhbF9tdXRzX2J5X2dyb3VwLnBkZiIsIHNlcCA9ICIiKSkKZ2dwbG90KHRvdGFsX211dGF0aW9uc190YWJsZV93bWV0YWRhdGEsIGFlcyh4ID0gYW5hbHlzaXMsIHkgPSB0b3RhbF9tdXRhdGlvbnNfc2NvcmUsIGZpbGwgPSBhbmFseXNpcykpICsgZ2VvbV9ib3hwbG90KG91dGxpZXIuc2hhcGUgPSBOQSkgKyAKICAgIGdlb21faml0dGVyKGFscGhhID0gMS84LCBhZXMoY29sb3VyID0gYW5hbHlzaXMpKSArIGd1aWRlcyhjb2xvdXIgPSBGQUxTRSwgZmlsbCA9IEZBTFNFKSArIHhsYWIoIkdyb3VwIikgKyB5bGFiKCJNdXRhdGlvbiBTY29yZSIpICsgCiAgICBnZ3RpdGxlKCJEaXN0cmlidXRpb24gb2YgTXV0YXRpb24gU2NvcmVzIHBlciBDZWxsIikgKyAKICAgIGdlb21fc2lnbmlmKGNvbXBhcmlzb25zID0gbGlzdChjKCJncm91cGVkX3BkIiwgImdyb3VwZWRfcHIiKSksIG1hcF9zaWduaWZfbGV2ZWw9VFJVRSwgeV9wb3NpdGlvbiA9IDQ1LCBhbm5vdGF0aW9ucyA9ICc8IDJlLTE2JykgKyAKICAgIGdlb21fc2lnbmlmKGNvbXBhcmlzb25zID0gbGlzdChjKCJncm91cGVkX3BkIiwgIm5haXZlIikpLCBtYXBfc2lnbmlmX2xldmVsPVRSVUUsIHlfcG9zaXRpb24gPSA1MCwgYW5ub3RhdGlvbnMgPSAnPCAyZS0xNicpICsgCiAgICBnZW9tX3NpZ25pZihjb21wYXJpc29ucyA9IGxpc3QoYygiZ3JvdXBlZF9wciIsICJuYWl2ZSIpKSwgbWFwX3NpZ25pZl9sZXZlbD1UUlVFLCB5X3Bvc2l0aW9uID0gNDAsIGFubm90YXRpb25zID0gJzAuMDIzJykgKwogICAgc2NhbGVfeF9kaXNjcmV0ZShsYWJlbHMgPSBjKCJuYWl2ZSIgPSAiVE4iLCAiZ3JvdXBlZF9wciIgPSAiUEVSIiwgImdyb3VwZWRfcGQiID0gIlBEIiksIGxpbWl0cz1jKCJuYWl2ZSIsImdyb3VwZWRfcHIiLCJncm91cGVkX3BkIikpCmRldi5vZmYoKQoKcGRmKGZpbGUgPSBwYXN0ZShkaXIsICJwbG90X291dC9OSTA4L05JMDhfdG90YWxfbXV0c19ieV9ncm91cF9yaWRnZXBsb3QucGRmIiwgc2VwID0gIiIpKQpnZ3Bsb3QodG90YWxfbXV0YXRpb25zX3RhYmxlX3dtZXRhZGF0YSwgYWVzKHggPSB0b3RhbF9tdXRhdGlvbnNfc2NvcmUsIHkgPSBhbmFseXNpcywgZmlsbCA9IGFuYWx5c2lzKSkgKyBnZW9tX2RlbnNpdHlfcmlkZ2VzKGFscGhhPS43KSArIAogICAgeWxhYigiR3JvdXAiKSArIHhsYWIoIk11dGF0aW9uIFNjb3JlIikgKyBzY2FsZV95X2Rpc2NyZXRlKGxhYmVscyA9IGMoIlBEIiwgIlBFUiIsICJUTiIpKSArIAogICAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKSArIHRoZW1lKGF4aXMudGl0bGUueSA9IGVsZW1lbnRfYmxhbmsoKSkKZGV2Lm9mZigpCgpwZGYoZmlsZSA9IHBhc3RlKGRpciwgInBsb3Rfb3V0L05JMDgvTkkwOF90b3RhbF9tdXRzX2J5X2dyb3VwX2FuZF9wYXRpZW50LnBkZiIsIHNlcCA9ICIiKSkKZ2dwbG90KHRvdGFsX211dGF0aW9uc190YWJsZV93bWV0YWRhdGEsIGFlcyh4ID0gYW5hbHlzaXMsIHkgPSB0b3RhbF9tdXRhdGlvbnNfc2NvcmUsIGNvbG9yID0gcGF0aWVudF9pZCkpICsKICAgIGdlb21fcG9pbnQocG9zaXRpb24gPSBwb3NpdGlvbl9qaXR0ZXJkb2RnZShkb2RnZS53aWR0aD0wLjkpKSArCiAgICBnZW9tX2JveHBsb3Qob3V0bGllci5jb2xvdXIgPSBOQSwgcG9zaXRpb24gPSBwb3NpdGlvbl9kb2RnZSh3aWR0aD0wLjkpKSArCiAgICB4bGFiKCJHcm91cCIpICsKICAgIHlsYWIoIk11dGF0aW9uIFNjb3JlIikgKwogICAgZ2d0aXRsZSgiRGlzdHJpYnV0aW9uIG9mIE11dGF0aW9uIFNjb3JlcyBwZXIgQ2VsbCIpICsKICAgIGd1aWRlcyhjb2xvdXI9Z3VpZGVfbGVnZW5kKHRpdGxlPSJSZXNwb25zZSBHcm91cCIpKSArCiAgICBndWlkZXMoZmlsbCA9IEZBTFNFLCBjb2xvciA9IEZBTFNFKSArCiAgICBzY2FsZV94X2Rpc2NyZXRlKGxhYmVscyA9IGMoIm5haXZlIiA9ICJUTiIsICJncm91cGVkX3ByIiA9ICJQRVIiLCAiZ3JvdXBlZF9wZCIgPSAiUEQiKSwgbGltaXRzPWMoIm5haXZlIiwiZ3JvdXBlZF9wciIsImdyb3VwZWRfcGQiKSkKZGV2Lm9mZigpCgpwZGYoZmlsZSA9IHBhc3RlKGRpciwgInBsb3Rfb3V0L05JMDgvTkkwOF90b3RhbF9tdXRzX2J5X3ByaW12c21ldC5wZGYiLCBzZXAgPSAiIikpCmdncGxvdChkYXRhPXN1YnNldCh0b3RhbF9tdXRhdGlvbnNfdGFibGVfd21ldGFkYXRhLCAhaXMubmEocHJpbWFyeV9vcl9tZXRhc3RhaWMpKSwgYWVzKHggPSBwcmltYXJ5X29yX21ldGFzdGFpYyAsIHkgPSB0b3RhbF9tdXRhdGlvbnNfc2NvcmUsIGNvbG9yID0gYW5hbHlzaXMpKSArCiAgICAgIGdlb21fcG9pbnQocG9zaXRpb24gPSBwb3NpdGlvbl9qaXR0ZXJkb2RnZShkb2RnZS53aWR0aD0wLjkpKSArCiAgICAgIGdlb21fYm94cGxvdChvdXRsaWVyLmNvbG91ciA9IE5BLCBwb3NpdGlvbiA9IHBvc2l0aW9uX2RvZGdlKHdpZHRoPTAuOSkpICsKICAgICAgeGxhYigiR3JvdXAiKSArCiAgICAgIHlsYWIoIk11dGF0aW9uIFNjb3JlIikgKwogICAgICBnZ3RpdGxlKCJEaXN0cmlidXRpb24gb2YgTXV0YXRpb24gU2NvcmVzIHBlciBDZWxsIikgKwogICAgICBndWlkZXMoY29sb3VyPWd1aWRlX2xlZ2VuZCh0aXRsZT0iUmVzcG9uc2UgR3JvdXAiKSkgKwogICAgICBndWlkZXMoZmlsbCA9IEZBTFNFKQpkZXYub2ZmKCkKYGBgCgo=